QUANTIFYING MODEL UNCERTAINTIES 
IN COMPLEX SYSTEMS 



JIARUI YANG AND JINQIAO DUAN 

Abstract. Uncertainties are abundant in complex systems. Appropriate mathematical models for 
these systems thus contain random effects or noises. The models are often in the form of stochas- 
tic differential equations, with some parameters to be determined by observations. The stochastic 
differential equations may be driven by Brownian motion, fractional Brownian motion, or Levy 
motion. 

After a brief overview of recent advances in estimating parameters in stochastic differential equa- 
tions, various numerical algorithms for computing parameters are implemented. The numerical 
simulation results are shown to be consistent with theoretical analysis. Moreover, for fractional 
Brownian motion and a- stable Levy motion, several algorithms are reviewed and implemented to 
numerically estimate the Hurst parameter H and characteristic exponent a. 



1. Introduction 

Since random fluctuations are common in the real world, mathematical models for complex sys- 
tems are often subject to uncertainties, such as fluctuating forces, uncertain parameters, or random 
boundary conditions Il89l 1551 1441 11211 11221 11251 . Uncertainties may also be caused by the lack 
of knowledge of some physical, chemical or biological mechanisms that are not well understood, 
and thus are not appropriately represented (or missed completely) in the mathematical models 

tm m [97i [mi w. 

Although these fluctuations and unrepresented mechanisms may be very small or very fast, their 
long-term impact on the system evolution may be delicate [[71 [551 SHI or even profound. This 
kind of delicate impact on the overall evolution of dynamical systems has been observed in, for 
example, stochastic bifurcation ll2~51[T7l 15511 . stochastic resonance |[59l , and noise-induced pattern 
formation ll44l[T4l . Thus taking stochastic effects into account is of central importance for mathe- 
matical modeling of complex systems under uncertainty. Stochastic differential equations (SDEs) 
are appropriate models for many of these systems f7l l27l[T08lll22L 

For example, the Langevin type models are stochastic differential equations describing various 
phenomena in physics, biology, and other fields. SDEs are used to model various price processes, 
exchange rates, and interest rates, among others, in finance. Noises in these SDEs may be mod- 
eled as a generalized time derivative of some distinguished stochastic processes, such as Brownian 
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motion (BM), Levy motion (LM) or fractional Brownian motion (fBM) [36]. Usually we choose 
different noise processes according to the statistical property of the observational data. For exam- 
ple, if the data has long-range dependence, we consider fractional Brownian motion rather than 
Brownian motion. If the data has considerable discrepancy with Gaussianity or normality, Levy 
motion may be an appropriate choice. In building these SDE models, some parameters appear, as 
we do not know certain quantities exactly. 

Based on the choice of noise processes, different mathematical techniques are needed in estimating 
the parameters in SDEs with Brownian motion, Levy motion, or fractional Brownain motion. 

In this article, we are interested in estimating and computing parameters contained in stochastic 
differential equations, so that we obtain computational models useful for investigating complex 
dynamics under uncertainty. We first review recent theoretical results in estimating parameters in 
SDEs, including statistical properties and convergence of various estimates. Then we develop and 
implement numerical algorithms in approximating these parameters. 

Theoretical results on parameter estimations for SDEs driven by Brownian motion are relatively 
well developed ((51|28l|48l|99l), and various numerical simulations for these parameter estimates 
([ffl[3l|99l[62]|) are implemented. So, in Section 2 below, we do not present such numerical results. 
Instead, we will concentrate on numerical algorithms for parameter estimations in SDEs driven by 
fractional Brownian motion and Levy motion in Section 3 and 4, respectively. 

This paper is organized as follows. In Section 2, we consider parameter estimation for SDEs with 
Brownian motion B t . We present a brief overview of some available techniques on estimating 
parameters in these stochastic differential equations with continuous-time or discrete-time obser- 
vations. In fact, we present results about how to estimate parameters in diffusion terms and drift 
terms, given continuous observations and discrete observations, respectively. 

In Section 3, we consider parameter estimation for SDEs driven by fractional Brownian motion 5f 
with Hurst parameter H. After discussing basic properties of fBM, we consider parameter estima- 
tion methods for Hurst parameter H from given fBM data. Then, we compare the convergence rate 
of each method by comparing estimates computed with hypothetic data. Unlike the case of SDEs 
with Brownian motion, there is no general formula for the estimate of the parameter in the drift 
(or diffusion) coefficient of a stochastic differential equation driven by fBM. We discuss different 
estimates associated with different models and discuss the statistical properties respectively. We 
develop and implement numerical simulation methods for these estimates. 

Finally, in Section 4, for stochastic differential equation with (non-Gaussian) a-stable Levy mo- 
tion Lf, we consider estimates and their numerical implementation for parameter a and other 
parameters in the drift or diffusion coefficients. 

2. Quantifying Uncertainties in SDEs Driven by Brownian motion 

In this section, we consider a scalar diffusion process X t e W l , < t < T satisfying the following 
stochastic differential equation 

( 1 ) dX t = fi(6, t, X t )dt + cr(#, t, X t )dB t , X = £ 

where B t is a m-dimensional Brownian motion, 6 e a compact subset of W and i?eSa compact 
subset of R. 9 are unknown parameters which are to be estimated on the basis of observations. Here 
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jx : 0x[OJ]x R d -> R d , the drift coefficient, and cr : S x [0, T]xW> ^ R dxm , the diffusion 
coefficient, are usually known functions but with unknown parameters 6 and &. 

Some remarks are in order here. 

• Under local Lipschitz and the sub-linear growth conditions on the coefficients /u and cr, 
there exists a unique strong solution of the above stochastic differential equation (see 117711 
or ll85l ) and this is an universal assumption for all results we discuss below. 

• The diffusion coefficient cr is almost surely determined by the process, i.e., it can be esti- 
mated without any error if observed continuously throughout a time interval (see RDl and 
1131). 

• The diffusion matrix defined by E(#, t, X t ) = cr(#, t, X t )cr{&, t, X t ) T plays an important role 
on parameter estimation problems. 

2.1. How to Estimate Parameters Given Continuous Observation. Since it is not easy to esti- 
mate parameters 6 and # at the same time, usually we simplify our model by assuming one of those 
parameters is known and then estimate the other. Moreover, instead of representing the results of 
all types of diffusion processes, we choose to present the conclusion of the most general one, such 
as, we prefer the nonhomogeneous case rather than the homogeneous one or the nonlinear one 
rather than the linear one. 

2.1.1. Parameter Estimation of Diffusion Terms with Continuous Observation. We assume that 
the unknown parameter 8 in the drift coefficient is known. Then our model can be simplified as 

(2) dX t = /d(t, X)dt + o-(&, t, X t )dB t , X Q = £ 
Remarks: 

• Different with the model (QQ), the drift coefficient /u(t,X) in model © is possibly unknown 
and maybe related to the whole past of process X instead of X,. In this case, our model can 
be easily extended to the non-Markovian case which is more general than case (1). 

• If fj, is depending on the unknown parameter # in model ©, we can also prove the local 
asymptotic mixed normality property for the maximum likelihood estimate (MLE) when 
H(t, X) = //(#, X t ) and cr(#, t, X t ) = cr(&, X,) (see I0S1 ). 

If the diffusion matrix £(#, t, X t ) is invertible, then define a family of contrasts by 

1 " 

(3) U n {§) = - V [log det fl v X r , ) + (X; i ) r E(#, f t _ v X,. , 

i=\ 

where 

XI = i(X,. - X ti ), 6'! = f t - t'U ,for\<i< n, 

and f- is an appropriate partition of [0,T]. However, this assumption does not always hold. So, we 
consider a more general class of contrasts of the form 

1 " 

(4) U\&) = - V tl v X fl ), XI), 
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where / should satisfy certain conditions to obtain the asymptotic property and consistency prop- 
erty for the estimate generated by these contrasts below (see II451I ). 

Let ■&„ be a minimum contrast estimate associated with U", i.e. § n satisfies the following equation 

U n (& n ) = min £/"(#). 

Under some smoothness assumptions on the coefficient /u and 6, empirical sampling measure as- 
sumption on the sample times t", and identifiability assumption on the law of the solution of ©, 
Genon-Catalot and Jacod IflTll have proved that the estimate § n has a local asymptotic mixed nor- 
mality property, i.e., V^($« _ #o) where # is the true value of the parameter converges in law to 
N(0, S). 

Remarks: 

• We do not include the drift coefficient /u in the contrast U"{d) because it is possibly un- 
known. Even if it is known, we still do not want it involved since it is a function of the 
whole past of X and thus is not observable. 

• If the diffusion matrix £ is invertible, it can be proven that the contrast of form © is optimal 
in the class of contrasts of type ©. 

2.1.2. Parameter Estimation of Drift Terms with Continuous Observations. We assume that the 
unknown parameter ■& in the diffusion coefficient is known. Then the model (OQ) can be simplified 

as 

(5) dX, = fj.(6, t, X t )dt + o-(t, X,)dB t , X = f . 

Since no good result for above general model exists, we introduce the result for the following non- 
homogeneous diffusion process instead. 

Consider a real valued diffusion process {X t , t > 0} satisfying the following stochastic differential 
equation: 

(6) dX t = fx(6, t, X,)dt + dB„ X = {, 

where the drift coefficient function /i is assumed to be nonanticipative. Denote the observation of 
the process by X^ := {X t , < t < T] and let P T e be the measure generated by the process X?. Then 
the Radon-Nicodym derivative (likelihood function) of P T g with respect to P^ where 6q is the true 
value of the parameter 6 is given by (see It80l0 

L T {6) := (dP T e ldP T d0 )(Xl) 

= exp{ f \fi(9, t, X t ) - fi(0 , t, X t )]dX t -if \j/(e, t, X t ) - (x 2 (8 , t, X t )]dt). 
Jo ^ Jo 

So we can get the Maximal Likelihood Estimate (MLE) defined by 

6 T := argsup ee ®L T (6). 

Then we can show that the MLE is strongly consistent, i.e., 6 T — » 9 P 6o - a.s. as T — » oo, and 
converge to a normal distribution (see Chapter 4 in lfT3ll for more details). 
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Remarks: 

• In [fl3l . Bishwal also proves that the MLE and a regular class of Bayes estimates (BE) are 
asymptotically equivalent. 

• By applying an increasing transformation as described in |fT1, 



r x du 

J cr(u) 



(7) Y t = g(X) 

- i) 

we can transform the diffusion process X, defined by 

dX t = /u(6,X t )dt + cr(X t )dB t 
into another diffusion process Y, defined by 

df, = p.{9, Y t )dt + dB t , 

where 



(8) my) = 



o-(g~ l (y)) 2 dy 



Then, we can get the MLE of process X, by calculating the MLE of process Y, according 
to what we learned in this section (see Q or for more details). 

2.2. How to Estimate Parameters given Discrete Observation. Given the practical difficulty in 
obtaining a complete continuous observation, we now discuss parameter estimations with discrete 
observation. 



2.2. 1 . Parameter Estimation of Drift Terms with Discrete Time. In this section, we assume that the 
unknown parameter # in the diffusion coefficient <x is known. Then the model (OQ) can be simplified 
as 

(9) dX t = fi(6, t, X t )dt + cr(t, X t )dB t , X = £. 

Ideally, when the transition densities p(s, x, t, y; 6) of X are known, we can use the log likelihood 
function 

n 

W) = log p(f M , X t ._ t , t u X t , ; 6), 
i=i 

to compute the LME 6 which is strongly consistent and asymptotically normally distributed, (see 
OH and JH, El and lfT09ll ). 

If the transition densities of X are unknown, instead of computing the log likelihood function /„(#), 
we would like to use approximate log-likelihood function which, under some regularity conditions 
(see |[56ll ), is given by 

Jo cr 2 (t, X t ) 2 Jo cr 2 (t, X t ) 
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to approximate the log-likelihood function based on continuous observations (see IU03H ). Then, 
using an ltd type approximation for the stochastic integral we can obtain 



^— ' <T (ti-\,X tj ,) 



1=1 



1 M (Q, ti-uX^) 
2j^ o-^ti-uXt.^) 



Thus, the maximizer of l n (6) provides an approximate maximum likelihood estimate (AMLE). In 
1992, Yoshida HI 3011 proved that the AMLE is weakly consistent and asymptotically normally dis- 
tributed when the diffusion is homogeneous and ergodic. In IfTBl . Bishwal got the similar result 
for the nonhomogeneous case with drift function t, X) = 9f(t, X t ) for some smooth functions 
f(t,x). Moreover, he measured the loss of information using several AMLEs according to different 
approximations to It(0). 



2.2.2. Parameter Estimation of Diffusion Terms (and/or Drift Terms) with Discrete Observation. 
In previous sections, we always assume one of those parameters is known and then estimate the 
other one. In this section, I want to include the situation when both 6 and # are unknown and how 
to estimate them based on the discrete observation of the diffusion process at the same time. 

Suppose we are considering the real valued diffusion process X t satisfying the following stochastic 
differential equation 

(10) dX t = ji(G, X,)dt + cr(#, X t )dB t . 

Denote the observation times by to = 0, t\, t 2 , . . . , t Nr , where Nt is the smallest integer such that 
t Nt+1 > T. In this section, we mainly consider three cases of estimating /3 = (6, jointly, = 9 
with § known and (3 = ■& with 6 known. In regular circumstances, the estimate /3 converges in 
probability to some ySand <T(fi-p) converges in law to A^(0, Qg) as T tends to infinity, where /3 
is the true value of the parameter. 

For simplicity, we set the law of the sampling intervals A„ =t„ — t„_i as 

(11) A = eA , 
where A has a given finite distribution and e is deterministic. 

Remark: We are not only studying the case when the sampling interval is fixed, i.e., Var[A ] = 0, 
but also the continuous observation case, i.e., 6 = and the random sampling case. 

Let h(y\, y , S,fi, e) denote a r-dimensional vector function which consists of r moment conditions 
of the discretized stochastic differential equation (flOl) (please see [f5TTl or ll54l for more details). 
Moreover, this function also satisfies 

E AM _ 1 [h(Y n ,Y n - 1 ,A n ,p,e)] = 0, 

where the expectation is taken with respect to the joint law of (A„, Y n , Y n _{). 
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By the Law of Large Numbers, E[h(Y n , Y n _i,A n ,/3, e)] may be estimated by the sample average 
defined by 

N T -1 

(12) m T {p) = Nj l Y,h{Y n ,Y n _ u K.p,e). 

n=l 

Then we can obtain an estimate j5 by minimizing a quadratic function 

(13) Q T (J3) = m T (j3)'W T m T (j3), 

where W T is a r x r positive definite weight matrix and this method is called Generalized Method 
of Moments (GMM). In [51], Hansen proved the strong consistency and asymptotic normality of 
GMM estimate, i.e. 

Vn0 - 0) -> N(0, V), 

when 6 = § and W T satisfied certain conditions. Mykland used this technique to obtain the closed 
form for the asymptotic bias but sacrificed the consistency of the estimate. 



3. Quantifying Uncertainties in SDEs Driven by Fractional Brownian Motion 

Colored noise, or noise with non-zero correlation in time, are common in physical, biological and 
engineering sciences. One candidate for modeling colored noise is fractional Brownian motion 

m. 

3.1. Fractional Brownian Motion. Fractional Brwonian motion (fBM) was introduced within a 
Hilbert space framework by Kolmogorov in 1940 in f73~l . where it was called Wiener Helix. It was 
further studied by Yaglom in H127II . The name fractional Brownian motion is due to Mandelbrot 
and Van Ness, who in 1968 provided in Il84l a stochastic integral representation of this process in 
terms of a standard Brownian motion. 

Definition 3.1 (Fractional Brownian motion D96|). Let H be a constant belonging to (0,1). A frac- 
tional Brownian motion (fBM) (B H (t)) t > of Hurst index H is a continuous and centered Gaussian 
process with covariance function 

By the above definition, we see that a standard fBM B H has the following properties: 

(1) 5^(0) = and E[B H (t)] = Ofor all t > 0. 

(2) B H has homogeneous increments, i.e., B H (t + s) — B H (s) has the same law of B H (t) for 
s,t>0. 

(3) B H is a Gaussian process and E[B H (t) 2 ] = t 2H , t < 0,for all H £ (0, 1). 

(4) B H has continuous trajectories. 

Using the method presented in ll23l |24| . we can simulate sample paths of fractional Brwonian 
motion with different Hurst parameters (see Figure CO). 

For H = 1/2, the fBM is then a standard Brownian motion. Hence, in this case the increments of 
the process are independent. On the contrary, for H ^ 1/2 the increments are not independent. 
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Figure 1. Three sample paths of fBM with Hurst parameter H = 0.25, 0.5, 0.9 

More precisely, by the definition of fBM, we know that the covariance between B H (t + h) - B H (t) 
and B H (s + h) - B H (s) with s + h <t and t- s = nh is 

p H (n)= 1 -h 2H [(n + ir + (n-l) 2H -2n 2H ]. 

In particular, we obtain that the two increments of the form B H (t+h)-B H (t) and B H (t+2h)- B H (t+h) 
are positively correlated for H > 1 /2, while they are negatively correlated for H < 1/2. In the first 
case the process presents an aggregation behavior and this property can be used in order to describe 
"cluster" phenomena (systems with memory and persistence). In the second case it can be used to 
model sequences with inter mitt ency and antipersistence. 

From the above description, we can get a general ideal that the Hurst parameter H plays an impor- 
tant role on how respective fBM behaves. So, it should be considered as an extra parameter when 
we estimate others in the coefficients of the SDE driven by fBM. 

Considering the further computation, we would like to introduce one more useful property of fBM. 

Definition 3.2 (Self-similarity). A stochastic process X= {X t , ( e I) is called 6-self-similar or 
satisfies the property of self-similarity if for every a > there exists b > such that 

Law(X at , t>0) = Law(a b X t , t > 0). 

Note that 'Law=Law' means that the two processes X at and a b X t have the same finite-dimensional 
distribution functions, i.e., for every choice to,...,t n in R, 

P(X ato < xq, . . . , X atn < x n ) = P(a b X tQ < xq, . . . , a b X tn < x n ). 

for every xq,...,x„ in R. 

Since the covariance function of the fBM is homogeneous of order 2H, we obtain that B H is a self- 
similar process with Hurst index H, i.e., for any constant a > 0, the processes B H {at) and a H B H (t) 
have the same distribution law. 

3.2. How to Estimate Hurst Parameter H. Let's start with the simplest case: 

dX t = dB H (t), i.e.,X, = B H (t) t > 0, 
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where {B H (t),t > 0} is a fBM with Hurst parameter H e (0, 1). Now, our question is how to 
estimate Hurst parameter H given observation data X , X\,..., X N . For a parameter estimation of 
Hurst parameter H, we need an extra ingredient, fractional Gaussian noise (fGn). 

Definition 3.3 (Fractional Gaussian noise). H110H 

Fractional Gaussian noise (fGn ) {Yj, i > 1 } is the increment of fractional Brownian motion, namely 

Y t = B H (i+ l)-B H (i), i>l. 

Remark: It is a mean zero, stationary Gaussian time series whose autocovariance function is given 
by 

p(h) = E(YjY i+h ) = + \) m - 2h 2H + \h- 1| 2// }, h > 0. 

An important point about p(h) is 

p(h) ~ H(2H - \)h 2H ~ 2 , as h -» oo, 

when H ± 1/2. Since p(h) = for h > 1 when H=l/2, the X,'s are white noise in this case. The 
Xi's, however, are positively correlated when \ < H < 1, and we say that they display long-range 
dependence (LRD) or long memory. 

From the expression of fGn, we know it is the same to estimate the Hurst parameter of fBM as 
to estimate the Hurst parameter of the respective fGn. Here, we introduce 4 different methods for 
measuring the Hurst parameter. Measurements are given on artificial data and the results of each 
method are compared in the end. However, the measurement techniques used in this paper can 
only be described briefly but references to fuller descriptions with mathematical details are given. 

3.2. 1. R/S Method. The R/S method is one of the oldest and best known techniques for estimating 
H. It is discussed in detail in (83 and EDI, p.83-87. 

For a time series {Y, : t = l,2,...,N} with partial sums given by Z(n) = YH=\ Yi and the sample 
variance given by 

S\n) = J—f]Yf- — !— Z(n) 2 , 
n - 1 j-t n(n - 1) 



then the R/S statistic, or the rescaled adjusted range, is given by: 

max(z(0 - -Z(n)) - min \Z(t) - -Z(ri)\ 

\<t<n \ n I l<t<n\ n I 



R 1 

77(") = 



5 S(n) 
For fractional Gaussian noise, 

E[R/S(n)] ~ C H n H , 
as n —* oo, where C H is another positive, finite constant not dependent on n. 

The procedure to estimate H is therefore as follows. For a time series of length N, subdivide the 
series into K blocks with each of size n = N/K. Then, for each lag n, compute RIS{k h ri), starting 
at points k t = iN/K + l,i = 1,2, . . . , K - 1. In this way, a number of estimates of R/S(n) are 
obtained for each value of n. For values of n approaching N, one gets fewer values, as few as 1 
when n>N- N/K. 
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Choosing logarithmically spaced values of n, plot \og[R/S(k h n)] versus logn and get, for each n, 
several points on the plot. This plot is sometimes called the pox plot for the R/S statistic. The 
parameter H can be estimated by fitting a line to the points in the pox plot. 

There are several disadvantages with this technique. Most notably, there are more estimates of the 
statistic for low values of n where the statistic is affected most heavily by short range correlation 
behavior. On the other hand, for high values of n there are too few points for a reliable estimate. 
The values between these high and low cut off points should be used to estimate H but, in practice, 
often it is the case that widely differing values of H can be found by this method depending on 
the high and low cut off points chosen. To modify the R/S statistic, we can use a weighted sum of 
autocovariance instead of the sample variance. Details can be found in [|82| . 

3.2.2. Aggregated Variance. Given a time series {Y t : t = 1,2, .. .,N}, divide this into blocks of 
length m and aggregate the series over each block. 

j km 

Y M(k):=- y Y { , k= 1,2, [N/m]. 

i=(k-\)m+\ 

We compute its sample variance, 

1 N/m 

y^Yim) _ Y(Y {m \k) - Yf. 

N/m j—* 

where 

Y _ £/=i 
N 

is the sample mean. The sample variance should be asymptotically proportional to m 2H ~ 2 for large 
N/m and m. Then, for successive values of m, the sample variance of the aggregated series is 
plotted versus m on a log-log plot. So we can get the estimate of H by computing the gradient of 
that log-log plot. However, jumps in the mean and slowly decaying trends can severely affect this 
statistic. One technique to combat this is to difference the aggregated variance and work instead 
with 

VarY^+V - VarYW. 

3.2.3. Variance of Residuals. This method is described in more detail in H101L Take the series 
{Y t : t = 1, 2, . . . , N} and divide it into blocks of length m. Within each block calculate partial 
sums: Z k (t) = Z;=(^7)m+i Yi, k = 1 . . .N/m, t = 1 . . .m. For each block make a least squares fit 
to a line at + but. Subtract this line from the samples in the block to obtain the residuals and then 
calculate their variance 

V k = - Y(Z k (t) -a k - b k tf. 

The variance of residuals is proportional to m 2H . For the proof in the Gaussian case, see H118L 
This variance of residuals is computed for each block, and the median (or average) is computed 
over the blocks. A log-log plot versus m should follow a straight line with a slope of 2H. 
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- Aggregated Variance Method 

- Peri odog ram Method 

- Variance of Residual Method 
Modified Aggregated Variance Method 
R/S 

actual value 



Figure 2. Numerical estimation of the Hurst parameter H of fBM: Actual value 
H = 0.65 



3.2.4. Periodogram. The periodogram is a frequency domain technique described in [|49ll . For a 
time series {Y t : t = 1,2, . . . ,N}, it is defined by 



1 

2nN 



I V 7 " 



where /I is the frequency. In the finite variance case, 1(A) is an estimate of the spectral density of 
Y, and a series with long-range dependence will have a spectral density proportional to |/l| 1_2// for 
frequencies close to the origin. Therefore, the log-log plot of the periodogram versus the frequency 
displays a straight line with a slope of 1-2H. 



3.2.5. Results on Simulated Data. In this subsection, we would like to use artificial data to check 
the robustness of above techniques and compare the result in the end. 

For each of the simulation methods chosen, traces have been generated. Each trace is 10,000 points 
of data. Hurst parameters of 0.65 and 0.95 have been chosen to represent a low and a high level 
of long-range dependence in data. From the Figure [2] and Figure [31 we can see that the Variance 
of Residual Method and R/S have the most accurate result. The Modified Aggregated Variance 
Method improved a little bit over the original one, but both of them still fluctuate too much. 



3.3. How to Estimate Parameters in SDEs Driven by fBM. After we discuss how to estimate 
the Hurst parameter of a series of artificial fBM data, now we want to concern how to estimate 
the parameters of the linear/nonlinear stochastic differential equation(s) driven by fBM. The coef- 
ficients in the stochastic differential equation could be deterministic or random, linear or nonlinear. 
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- Aggregated Variance Method 

- Periodogram Method 

- Variance of Residual Method 
Modified Aggregated Variance Method 
R/S 

actual value 



1 



1.1 1.2 1.3 



1-5 1-6 



1.8 1.S 



Figure 3. Numerical estimation of the Hurst parameter H of fBM: Actual value 
H = 0.95 

No general results are available. So some specific statistical results will be discussed below ac- 
cording to what kind of specified models we deal with. 



3.3.1. Preparation. The main difficulty in dealing with a fBm is that it is not a semimartingale 
when H + \ and hence the results from the classical stochastic integration theory for semimartin- 
gales can not be applied. So, we would like to introduce the following integral transformation 
which can transform fBM to martingale firstly and it will be a key point in our development below. 
For < s < t < T, denote 



(14) 
(15) 

(16) 
(17) 



k H (t, s) 



-1 Jl/2)-H, 



,(1/2)-// 



2HY(3 /2 - H)Y(H + 1/2), 



h ,-i,2-2* , 2HT(3-2H)T(H+l/2) 
w, = A H t ; A H = — 



Mf = [ k H {t,s)dB^. 
Jo 



r(3/2 - H) 



Then the process M H is a Gaussian martingale (see f78l and 11921 ). called the fundamental martin- 
gale with variance function w H . 



3.3.2. Parameter Estimation for a Fractional Langevin Equation. Suppose {X t , t> 0} satisfies the 
following stochastic differential equation 



X t = 6 f X s ds + o-B?; 0<t<T, 
Jo 
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where 8 and <x 2 are unknown constant parameters, Bf is a fBM with Hurst parameter H e [1/2, 1]. 
Denote the process Z=(Z t , t e [0, T]) by 

(18) Z t = f k H (t,s)dX s . 

Jo 

Then we can prove that process Z is a semimartingale associated to X with following decomposi- 
tion (see 11691 ) 

(19) Z, = 6 f Q(s)dwf + <tM? , 

Jo 

where 

(20) Q(t) = T^H f k H (t,s)X(s)ds, 

dw" Jo 

and Mfh the Gaussian martingale defined by (17). From the representation (fT9l ), we know the 
quadratic variation of Z on the interval [0, t] is nothing but 

(Z) t = crw t , a.s. 

Hence the parameter cr 2 can be obtained by 

[wf r 1 lim J] K, - z ';f = a - s - 

i 

where t" is an appropriate partition of [0,t] such that sup ; \ f i+l - t"\ — » as n -* oo. So, the variance 
parameter can be computed with probability 1 on any finite time interval. 



As for the parameter 6, by applying the Girsanov type formula for fBM which is proved in 11691 , 
we can define the following maximum likelihood estimate of 6 based on the observation on the 
interval [0, t] by 

(21) ° T = {j Q2(s)dw ") X Q(s)dZs ' 

where processes Q, Z and are defined by (|20l ), (fl"8l and (16), respectively. For this estimate, 
strong consistency is proven and explicit formulas for the asymptotic bias and mean square error 
are derived by Kleptsyna and Le Breton 11701 . 
Remarks: 

• When H = 1/2, since Q = Z = X and dojj 1 = ds, the formula (I2TI) reduces to the result of 
Il80l for an usual Ornstein-Uhlenbeck process. 

• For an arbitrary He [1/2, 1], we could derive the following alternative expression for 6 T : 

\~ l ( i r-T 

l T s 2H ~ l dZ s -t\ 



*.{2j[ <*WHf} 



Example 3.4. Consider a special Ornstein-Uhlenbeck model 

dX t = 6X t dt + 2dBf . 
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Figure 4. Numerical estimation of drift parameter 9 in dX t = 6X,dt + 2dBf with 
Hurst parameter H = 0.75: Actual value 9 = 1 

Then, according to the above approximation scheme, we can numerically estimate 6=1 and the 
results are shown in Figure\4\ 

3.3.3. Parameter Estimation in Linear Deterministic Regression. Suppose X t satisfies the follow- 
ing stochastic differential equation 



X t = 6 f A(s)ds + f C(s)dB^, < t < T, 
Jo Jo 



where A and C are deterministic measurable functions on [0,T], Bf is a fBM with Hurst parameter 
H 6 [1/2, 1]. 



Let q t be defined by 



d f A 
<lt=—z k H (t,s)—(s)ds, 
dwy Jo C 



where wf and k H (t, s) are defined by (16) and CHI). Then, from Theorem 3 in ll69ll . we obtain the 
maximum likelihood estimate of 6 defined by 



6t = I J ^dw"j q t dZ t , 



where Z t is defined by ([18 



Remark: This result can be extended to an arbitrary H in (0,1) (see If78l ) and Q T is also the best 
linear unbiased estimate of 6. 

Example 3.5. Consider a special Linear Deterministic Regression 

dX t = -Bdt + tdBf. 
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Figure 5. Numerical estimation of drift parameter 9 in a Linear Deterministic Re- 
gression dX t = -9dt + tdBf with Hurst parameter H = 0.75: Actual value 6=1 



Then, using the above estimate, we can do numerical simulation with result shown in Figure\5\ 

3.3.4. Parameter Estimation in Linear Random Regression. Let us consider a stochastic differen- 
tial equation 

X(f) = [A(t,X(t)) + 6C(t,X(t))]dt + cr(t)dBf, t > 0, 

where B = {5f , t > 0} is a fractional Brownian motion with Hurst parameter H and o~(t) is a 
positive nonvanishing function on [0, oo). According to H105L the Maximum Likelihood Estimate 
9 T of 6 is given by 

_ £j 2 (t)dZ t + g WhWwf 
£j 2 2 {t)dw? 

where the processes Z t , J\,J 2 are defined by 

Z, = — — — dX s ,t>0, 
Jo o-(s) 

J w) = ~T~h k H<,t,s) — — ds,J 2 (t) = —— k H (t,s) — — ds, 

dw? Jo o-(s) dwy Jo o-(s) 

and wf, k H (t, s) are defined by (16) and (14). Also in the same paper, they proved that 6 T is strongly 
consistent for the true value 9. 

Example 3.6. Consider a special Linear Random Regression 

dX t = (t + 9X t )dt + tdBf. 
A numerical estimation of the parameter 9 is shown in Figure® 
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Figure 6. Numerical estimation of drift parameter 6 in a Linear Random Regression 
dX t = (t + 6X t )dt + tdB" with Hurst parameter H = 0.75: Actual value 6=1 

4. Parameter Estimation for SDE Driven by cc-Stable Levy Motion 

Brownian motion, as a Gaussian process, has been widely used to model fluctuations in engi- 
neering and science. For a particle in Brownian motion, its sample paths are continuous in time 
almost surely (i.e., no jumps), its mean square displacement increases linearly in time (i.e., normal 
diffusion), and its probability density function decays exponentially in space (i.e., light tail or ex- 
ponential relaxation) [95]. However some complex phenomena involve non-Gaussian fluctuations, 
with properties such as anomalous diffusion (mean square displacement is a nonlinear power law 
of time) |[T5l and heavy tail (non-exponential relaxation) H129H . For instance, it has been argued 
that diffusion in a case of geophysical turbulence H114H is anomalous. Loosely speaking, the diffu- 
sion process consists of a series of "pauses", when the particle is trapped by a coherent structure, 
and "flights" or "jumps" or other extreme events, when the particle moves in a jet flow. Moreover, 
anomalous electrical transport properties have been observed in some amorphous materials such 
as insulators, semiconductors and polymers, where transient current is asymptotically a power law 
function of time H112[[53l . Finally, some paleoclimatic data [29] indicates heavy tail distributions 
and some DNA data H114H shows long range power law decay for spatial correlation. Levy mo- 
tions are thought to be appropriate models for non-Gaussian processes with jumps HI 1 IB . Here we 
consider a special non-Gaussian process, the a-stable Levy motion, which arise in many complex 
systems H1261 . 

4.1. or-Stable Levy Motion. There are several reasons for using a stable distribution to model a 
fluctuation process in a dynamical system. Firstly, there are theoretical reasons for expecting a non- 
Gaussian stable model, e.g. hitting times for a Brownian motion yielding a Levy distribution, and 
reflection off a rotating mirror yielding a Cauchy distribution. The second reason is the Generalized 
Central Limit Theorem which states that the only possible non-trivial limit of normalized sums of 
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i.i.d. terms is stable. The third argument for modeling with stable distributions is empirical: Many 
large data sets exhibit heavy tails and skewness. In this section, we consider one-dimensional 
a-stable distributions defined as follows. 

Definition 4.1. fll64ll . Chapter 2.4) The Characteristic Function ip(u) of an a-stable random vari- 
able is given by 

(22) cp(u) = exp((-o- a )\u\ a {l - i/3sgn(u) tan(o7r/2)} + ifiu) 
where a e (0, 1) U (1,2), /3 e [-1, l],<x e R + , p e R, and by 

2 

(23) ip(u) = exp(-<x|w|{l + i/3-sgn(u)\og(\u\)} + ipu) 

n 

when a = 1, it gives a very well-known symmetric Cauchy distribution and 

1 , 

(24) (p(u) = exp(--cr|w| z + ipu), 
when a = 2, it gives the well-known Gaussian distribution. 

For the random variable X distributed according to the rule described above we use the notation 
X ~ S a (cr,/3,p). Especially when p = f5 = 0, i.e., X is a symmetric a-stable random variable, we 
will denote it as X ~ S aS . 

Also, from above definition, it is easy to see that the full stable class is characterized by four 
parameters, usually designated a,fi, cr, and fi. The shift parameter /i simply shifts the distribution to 
the left or right. The scale parameter cr compresses or extends the distribution about /j. in proportion 
to cr which means, if the variable x has the stable distribution X ~ S Q ,(cr, J S,/i), the transformed 
variable z = (x - fj)/cr will have the same shaped distribution, but with location parameter and 
scale parameter 1. The two remaining parameters completely determine the distribution's shape. 
The characteristic exponent a lies in the range (0, 2] and determines the rate at which the tails of 
the distribution taper off. When a = 2, a normal distribution results. When a < 2, the variance is 
infinite. When a > 1, the mean of the distribution exists and is equal to [i. However, when a < 1, 
the tails are so heavy that even the mean does not exist. The fourth parameter J3 determines the 
skewness of the distribution and lies in the range [-1,1]. 
Now let us introduce cr- stable Levy motions. 

Definition 4.2. (a-stable Levy motion 11641 ) 

A stochastic process {X(t) : t > 0} is called the (standard) a-stable Levy motion if 

(1) X(0)=0a.s.; 

(2) {X(t) : t > 0} has independent increments; 

(3) X(t)-X(s)~ S a ((t- s) l/a ,j3,0)forany0< s < t < oo. 

So, from the third condition, we can simulate all a-stable Levy motion if we know how to simulate 
X ~ S a (cr,fi, 0). Especially, it is enough to simulate X ~ S a (cr, 0, 0) if we want to get the trajecto- 
ries of symmetric a-stable Levy motions. 

We recall an important property of a-Stable random variables giving us the following result: It is 
enough to know how to simulate X ~ S a (l, 0, 0) in order to get any X ~ S a (a, 0, 0), Vcr e R + . 
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Figure 7. Three sample paths of symmetric a- stable Levy motion with a = 
0.4, 1.2, 1.9, respectively 

Proposition 4.3. If we have Xi,X 2 ~ S a (cr,8,p) and A, B are real positive constants and C is a 
real constant, then 

AX i +BX 2 + C ~ S a (o-(A a + B a ) 1/a ,B,p(A a + B a ) 1/a + C) 

Proposition 4.4. Let X ~ S a (o;B,0), with < a < 2, Then E\X\ P < oo far any < p < a, 

E\X\ P = oo far any p > a. 

Figure [7J shows sample paths of the a-stable Levy motion with different a. 

As we can see in Figure[71 the bigger the parameter a is, the more the path looks like Brownian mo- 
tion. Generally speaking, when we deal with concrete data, we have to choose a-stable processes 
very carefully to get the best estimation. We now discuss how to estimate a. 

4.2. How to Estimate the Characteristic Exponent a. Five different methods about how to esti- 
mate the characteristic exponent a of or-stable distribution are considered: Characteristic Function 
Method(CFM), Quantile Method, Maximum Likelihood Method, Extreme Value Method and Mo- 
ment Method. As in the last section, measurements are given on artificial data and the results of 
each method are compared in the end of this section. 

4.2.1. Characteristic Function Method. Since a-stable distributions are uniquely determined by 
their Characteristic Function (CF), it is natural to consider how to estimate parameter by studying 
their CF. Press H106II introduced a parameter estimation method based on CF, which gets estima- 
tions of parameters by minimizing differences between values of sample CF and the real ones. But 
this method is only applicable to standard distributions. 

Another method which uses the linearity of logarithm of CF was developed by Koutrouvelis [1741 
and it can be applied to general a-stable cases. This method is denoted as Kou-CFM. The idea is 
as follows: On the one hand, taking the logarithm of real part of CF gives 

\n[-Re((f(u))] = arln|w| + aTncr. 

On the other hand, the sample characteristic function of <p(Q) is given by ®(#) = (2f=i el6>k ) where 
Vfe's are n independent observations. In [74], a regression technique is applied to gain estimates 
for all parameters of a observed a stable distribution. In [1721 . Kogon improved this method by 
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replacing a linear regression fit by a linear least square fit which gave a more accurate estimate and 
its computational complexity became lower. 

4.2.2. Quantile Method. Quantiles are points taken at regular intervals from the cumulative distri- 
bution function of a random variable. Suppose we have n independent symmetric a-stable random 
variables with the stable distribution S a (cr,f3,(j.), whose parameters are to be estimated. Let x p be 
the p-th quantile, so that S a (x p ; ct,/3,ij.) = p. Let x p be the corresponding sample quantile, then x p 
is a consistent estimate of x p . 

In 1971, Fama and Roll [41] discovered that, for some large p (for example, p=0.95), 

= hZpZL = (0,827) jg I V 

2cr X J2 ~ *0.28 

is an estimate of the p-quantile of the standardized symmetric stable distribution with exponent a. 
According to this, they proposed a estimate (QM) for symmetric a-stable distributions. However, 
the serious disadvantage of this method is that its estimations are asymptotically biased. 
Later on, McCulloch Il871l improved and extended this result to general or-stable distributions, de- 
noted as McCulloch-QM. Firstly, he defined 

V a = (-X-0.95 _ X-ox)5)/(X-ojs — X_o.25) 

Vfi = (X-0.95 + *-0.05 _ 2jCo.5)/(A:-0.95 _ *-0.05) 

and let v a and vp be the corresponding sample value: 

Va = (*-0.95 _ *-0.05)/C*-0.75 _ *-0.25) 

Vp — (*-0.95 + -£-0.05 _ 2xo.5)/(x_o.95 _ #-0.05) 

which are the consistent estimates of the index v a and vp. Then, he illustrated that estimates of 
a can be expressed by a function of v a and Vp. Compared with QM, McCulloch-QM could get 
consistent and unbiased estimations for the general a-stable distribution, and extend the estimation 
range of parameter a to 0.6 < a < 2. Despite its computational simplicity, this method has a 
number of drawbacks, such as, there are no analytical expressions for the value of the fraction, and 
the evaluation of the tables implies that it is highly dependent on the value of a in a nonlinear way. 
This technique does not provide any closed-form solutions. 



4.2.3. Extreme Value Method. In 1996, based on asymptotic extreme value theory, order statistics 
and fractional lower order moments, Tsihrintzis and Nikias II 1 1911 proposed a new estimate which 
can be computed fast for symmetric a stable distribution from a set of i.i.d. observations. Five 
years later, Kuruoglu ll76l extended it to the general a stable distributions. The general idea of this 
method is as follows. Given a data series {Xj : i = 1, 2, . . . , N}, divide this into L nonoverlapping 
blocks of length K such that K = N/L. Then the logarithms of the maximum and minimum 
samples of each segment are computed as follows 

Yi = log(max{X lK - K+i \i = 1,2, . . . ,K}), 
Y, = log(-mm{X M+/ |/ = l,2,...,*}). 
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The sample means and variances of Y/ and Y t are calculated as 



1=1 1=1 



/=i z=i 
Finally, an estimate for a is given by sample variance as follows 

^ 2Vf>U £ 

Even though the accuracy and computational complexity decrease, there is now a closed form for 
the block size which means we have to look-up table to determine the segment size K. 

4.2.4. Moment Estimation Method. Another way to estimate parameters of the general a-stable 
distribution is the Logarithmic Moments Methods which was also introduced by Kuruoglu 117611 . 
The advantage of this method relative to the Fractional Lower Order Method is that it does not 
require the inversion of a sine function or the choice of a moment exponent p. The main feature is 
that the estimate of a can be expressed by a function of the second-order moment of the skewed 
process, i.e. 

-1/2 

& = { — 



J2 



where \j/ x = ^ and, for any X ~ S a (cr,fi, 0), L 2 is defined as follows 

1 1 

. 2 a 2 J a 



L 2 = £[(log \X\ - E[\og \X\]Y] = «Ai | - + — - — • 



4.2.5. Results on Simulated Data. In this subsection, we would like to use artificial data to check 
the robustness of the above techniques and compare the results. 



For each of the simulation methods chosen, estimates of a have been generated respectively and 
each trace is 1,000 points of data. Characteristic exponents of 0.95 and 1.70 have been chosen to 
represent a low and a high level of the rate at which the tails of the distribution taper off. 
From the Figures [8] and |9l we can see that the Characteristic Function Method and the Moment 
Estimate Method have the most accurate result. The Quantile Method behaved a little better than 
Extreme Value Method, but both of them still fluctuate too much when a is small. As to the 
convergence, we can see that all the methods get closer and closer to the real value when the points 
of data increase except for the Extreme Value Method. 

4.3. How to Estimate Parameters in SDEs Driven by Levy Motion. After we discussed how to 
estimate the characteristic exponent of a-stable Levy motions, now we consider how to estimate 
the parameters in stochastic differential equations driven by general Levy motion. Just as what 
we discussed about fBM, no general results about the parameter estimation for general cases are 
available at this time. Some special results will be listed below for different equations. 
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Figure 8. Numerical estimation of the characteristic exponent a in the a-stable 
Levy motion Lf by 4 different methods: Actual value a = 0.95 
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Figure 9. Numerical estimation of the characteristic exponent a in the a-stable 
Levy motion L" by 4 different methods: Actual value a = 1.70 



We consider parameter estimation of the Levy-driven stationary Ornstein-Uhlenbeck process. Re- 
cently, Brockdwell, Davis and Yang [16J studied parameter estimation problems for Levy-driven 
Langevin equation (whose solution is called an Ornstein-Uhlenceck process) based on observa- 
tions made at uniformly and closely-spaced times. The idea is to obtain a highly efficient estimate 
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of the Levy-driven Ornstein-Uhlenceck process coefficient by estimating the corresponding coef- 
ficient of the sampled process. The main feature is discussed below. 

Consider a stochastic differential equation driven by the Levy motion {L(t), t > 0} 

dY(t) = -9Y(t)dt + crdL(t). 
When L(t) is Brownian motion, the solution of above equation can be expressed as 



(25) Y(t) = e~ et Y(0) + cr f e~ 8(t ~ u) dL(t). 

Jo 

For any second-order driving Levy motion, the process {Y(t)} can be defined in the same way, and 
if [Lit)} is non-decreasing, {Y(t)} can also be defined pathwise as a Riemann-Stieltjes integral by 
(|25l ). For the convenience of the simulation, we rewrite solution as follows 



(26) Y(t) = e~ &{, ~ s) Y(s) + cr £ e~ e(t - u) dL(u), for all t > s>0. 

Now we collect all information corresponding to the sampled process in order to get the estimate. 
Set t = nh and s = (n - \)h in equation (|26l . Then the sampled process {Y^\ n = 0, 1,2, . . .} (or 
the discrete-time AR(1) process) satisfies 

where 

(p = e~ eh , andZ n = cr e- e(nh - u) dL(u). 

J(n-\)h 

Then, using the highly efficient Davis-McCormick estimate of <p, namely 

vQi) 

<P N = mlfl i<n<N-^, 



we can get the estimate of 6 and cr as follows 



.(2) zv n VrvW v (h \ 2 
a N - —jrZs- 1 N ' " 



i=0 

Example 4.5. Consider a Levy-driven Ornstein-Uhlenbeck process satisfying the following SDE 

(27) dX t = -X t dt + o-dL* . 

A numerical estimation of the diffusion parameter cr is shown in FigureUOl 
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Figure 10. Numerical estimation of the diffusion parameter <x in Levy-driven 
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£T = 2 



[4] S. Albeverrio, B. Rudiger and J. L. Wu (2000), Invariant measures and symmetry property of Levy type 

operators, Potential Analysis 13, 147-168. 
[5] S. Alizadeh, M. W. Brandt and F. X. Diebold (2002), Range-based estimation of stochastic volatility models, 

The Journal of Finance 57(3), 1047-1091. 
[6] D. Applebaum (2009), Levy Processes and Stochastic Calculus, 2nd edition, Cambridge University Press, UK. 
[7] L. Arnold (1998), Random Dynamical Systems, Springer, New York. 

[8] O. E. Barndorff-Nielsen, T. Mikosch and S. I. Resnick (Eds.) (2001), Levy Processes: Theory and Applications, 
Birkhauser, Boston. 

[9] C. Bender (2003), An Ito formula for generalized functionals of a Fractional Brownian motion with arbitrary 

Hurst parameter, Stock. Proc. Appl. 104, 81-106. 
[10] J. Beran (1994), Statistics for Long-Memory Processes, Chapman and Hall. 
[11] J. Bertoin (1998), Levy Processes, Cambridge University Press, Cambridge, U.K.. 
[12] P. Billingsley (1961), Statistical Inference for Markov Processes, Chicago University Press, Chicago. 
[13] J. P. N. Bishwal (2007), Parameter Estimation in Stochastic Differential Equations, Springer, New York. 
[14] D. Blomker and S. Maier-Paape (2003), Pattern formation below criticality forced by noise, Z. Angew. Math. 

Phys. 54(1), 1-25. 

[15] J. P. Bouchaud and A. Georges (1990), Anomalous diffusion in disordered media: Statistic mechanics, models 

and physical applications, Phys. Repts 195, 127-293. 
[16] P. J. Brockwell, R. A. Davis, and Y. Yang (2007), Estimation for nonnegative Levy-driven Ornstein-Uhlenbeck 

processes, J. Appl. Probab. 44(4), 977-989. 
[17] T. Caraballo, J. Langa and J. C. Robinson (2001), A stochastic pitchfork bifurcation in a reaction-diffusion 

equation, R. Soc. Lond. Proc. Ser. A 457, 2041-2061. 
[18] B. Chen (2009), Stochastic dynamics of water vapor in the climate system, Ph.D. Thesis, Illinois Institute of 

Technology, Chicago, USA. 

[19] B. Chen and J. Duan (2009), Stochastic quantification of missing mechanisms in dynamical systems, In "Recent 
Development in Stochastic Dynamics and Stochastic Analysis", Interdisciplinary Math, Sci. 8, 67-76. 



24 JIARUI YANG AND JINQIAO DUAN 

[20] A. Chronopoulou and F. Viens (2009), Hurst index estimation for self-similar processes with long-memory. In 
"Recent Development in Stochastic Dynamics and Stochastic Analysis", J. Duan, S. Luo and C. Wang (Eds.), 
World Scientific. 

[21] J. M. Corcuera, D. Nualart, and J. H. C. Woerner (2006), Power variation of some integral fractional processes, 
Bernoulli, 12,713-735. 

[22] J. M. Corcuera, D. Nualart and J. H. C. Woerner (2007), A functional central limit theorem for the realized 

power variation of integrated stable process, Stochastic Analysis and Applications 25, 169-186. 
[23] J. Coeurjolly (2001), Estimating the parameters of the Fractional Brwonian motion by discrete variations of its 

sample paths, Statistical Inference for Stochastic Processes 4, 199-227. 
[24] J. Coeurjolly (2000): Simulation and identification of the Fractional Brwonian motion: a bibliographical and 

comparative study, Journal of Statistical Software, American Statistical Association 5(07). 
[25] H. Crauel and F. Flandoli (1998), Additive noise destroys a pitchfork bifurcation, Journal of Dynamics and 

Differential Equations 10, 259-274. 
[26] D. Dacunha-Castelle adn D. Florens-Zmirou (1986), Estimation of the coefficients of a diffusion from discrete 

observations, 19, 263-284. 

[27] G. Da Prato and J. Zabczyk (1992), Stochastic Equations in Infinite Dimensions, Cambridge University Press. 

[28] M. Davis (2001), Pricing weather derivatives by marginal value, Quantitative Finance 1(3), 305-308. 

[29] P. D. Ditlevsen (1999), Observation of a-stable noise induced millennial climate changes from an ice record, 

Geophys. Res. Lett. 26, 1441-1444. 
[30] J. L. Doob (1953), Stochastic Processes, John Wiley, New York. 

[31] A. Du and J. Duan (2009), A stochastic approach for parameterizing unresolved scales in a system with memory, 

Journal of Algorithms & Computational Technology 3, 393-405. 
[32] J. Duan (2009), Stochastic modeling of unresolved scales in complex systems, Frontiers of Math, in China, 4, 

425-436. 

[33] J. Duan (2009), Predictability in spatially extended systems with model uncertainty I & II, Engineering Simu- 
lation!, 17-32 & 3 21-35. 

[34] J. Duan (2009), Predictability in nonlinear dynamical systems with model uncertainty, Stochastic Physics and 
Climate Modeling, T. N. Palmer and P. Williams (eds.), Cambridge Univ. Press, pp. 105-132. 

[35] J. Duan, X. Kan and B. Schmalfuss (2009), Canonical sample spaces for stochastic dynamical systems, In 
"Perspectives in Mathematical Sciences" , Interdisciplinary Math. Sci. 9, 53-70. 

[36] J. Duan, C. Li and and X. Wang (2009), Modeling colored noise by Fractional Brownian motion, Interdisci- 
plinary Math. Sci. 8, 119-130. 

[37] J. Duan and B. Nadiga (2007), Stochastic parameterization of large Eddy simulation of geophysical flows, Proc. 
American Math. Soc. 135, 1187-1196. 

[38] G. Dohnal (1987), On estimating the diffusion coefficient, J. Appl. Prob. 24, 105-1 14. 

[39] O. Elerian, S. Chib and N. Shephard (2001), Likelihood inference for discretely observed non-linear diffusions, 

Econometrica 69(4), 959-993. 
[40] B. Eraker (2001), MCMC analysis of diffusion models with application to finance, Journal of Business and 

Economic Statistics 19(2), 177-191. 
[41] E. F. Fama, R. Roll (1971), Parameter estimates for symmetric stable distribution Journal of the American 

Statistical Association, 66, 331-338. 
[42] D. Florens-Zmirou (1989), Approximate discrete-time schemes for statistics of diffusion processes, Statistics 

20, 547-557. 

[43] C. W. Gardiner (1985), Handbook of Stochastic Methods, Second Ed., Springer, New York. 

[44] J. Garcia-Ojalvo and J. M. Sancho (1999), Noise in Spatially Extended Systems, Springer- Verlag, 1999. 

[45] V. Genon-Catalot and J. Jacod (1993), On the estimation of the diffusion coefficient for multi-dimensional 
diffusion processes, Annates de VI. H. P., section B, tome 29, 1993. 

[46] V. Genon-Catalot and J. Jacod (1993), On the estimation of the diffusion coefficient for multidimensional diffu- 
sion processes, Ann. Inst. Henri Poincare, Probability's et Statistiques. 29, 119-151. 

[47] V. Genon-Catalot and J. Jacod (1994), On the estimation of the diffusion coefficient for diffusion processes, J. 
Statist. 21, 193-221. 



QUANTIFYING UNCERTAINTIES IN COMPLEX SYSTEMS 



25 



[48] V. Genon-Catalot, T. Jeantheau and C. Laredo (1999), Parameter estimation for discretely observed stochastic 

volatility models, Bernoulli 5(5), 855-872. 
[49] J. Geweke and S. Porter-Hudak (1983), The estimation and application of long memory time series models, 

Time Ser. Anal. 4, 221-238. 

[50] P. Hanggi and P. Jung (1995), Colored noise in dynamical systems, Advances in Chem. Phys. 89, 239-326. 
[51] L. P. Hansen (1982), Large sample properties of generalized method of moments estimators, Econometrica 63, 
767-804. 

[52] C. Hein, P. Imkeller and I. Pavlyukevich (2009), Limit theorems for p-variations of solutions of SDEs driven by 
additive stable Levy noise and model selection for paleo-climatic data, In "Recent Development in Stochastic 
Dynamics and Stochastic Analysis", J. Duan, S. Luo and C. Wang (Eds.), Interdisciplinary Math. Sci. 8. 

[53] M. P. Herrchen (2001), Stochastic Modeling of Dispersive Diffusion by Non-Gaussian Noise, Doctorial Thesis, 
Swiss Federal Inst, of Tech., Zurich. 

[54] C. C. Heyde (1997), Quasi-Likelihood and its Application: A General Approach to Optimal Parameter Estima- 
tion. Springer, New York. 

[55] W. Horsthemke and R. Lefever (1984), Noise-Induced Transitions, Springer- Verlag, Berlin. 

[56] J. E. Hutton and P. I. Nelson (1986), Quasi-likelihood estimation for semimartingales, Stochastic Processes and 
their Applications 22, 245-257. 

[57] I. A. Ibragimov, R. Z. Has'minskii (1981), Statistical Estimation-Asymptotic Theory. Springer- Verlag. 

[58] N. Ikeda and S. Watanabe (1989), Stochastic Differential Equations and Diffusion Processes, North-Holland 
Publishing Company, Amsterdam. 

[59] P. Imkeller and I. Pavlyukevich (2002), Model reduction and stochastic resonance, Stochastics and Dynamics 
2(4), 463-506. 

[60] P. Imkeller and I. Pavlyukevich (2006), First exit time of SDEs driven by stable Levy processes, Stoch. Proc. 
Appl. 116,611-642. 

[61] P. Imkeller, I. Pavlyukevich and T. Wetzel (2009), First exit times for Levy-driven diffusions with exponentially 
light jumps, Annals of Probability 37(2), 530C564. 

[62] J. Nicolau (2004), Introduction to the Estimation of Stochastic Differential Equations Based on Discrete Obser- 
vations, Stochastic Finance 2004 (Autumn School and International Conference). 

[63] J. Jacod (2006), Parametric inference for discretely observed non-ergodic diffusions, Bernoulli 12(3), 383-401. 

[64] A. Janicki and A. Weron (1994), Simulation and Chaotic Behavior of a-Stable Stochastic Processes, Marcel 
Dekker, Inc.. 

[65] W. Just, H. Kantz, C. Rodenbeck and M. Helm (2001), Stochastic modeling: replacing fast degrees of freedom 
by noise, /. Phys. A: Math. Gen. 34, 3199-3213. 

[66] I. Karatzas and S. E. Shreve (1991), Brownian Motion and Stochastic Calculus 2nd edition, Springer. 

[67] M. Kessler (2000), Simple and explicit estimating functions for a discretely observed diffusion process, Scan- 
dinavian Journal of Statistics 27(1), 65-82. 

[68] V. Krishnan (2005), Nonlinear Filtering and Smoothing: An Introduction to Martingales, Stochastic Integrals 
and Estimation, Dover Publications, Inc., New York. 

[69] M. L. Kleptsyna, A. Le Breton and M. C. Roubaud (2000), Parameter estimation and optimal filtering for 
fractional type stochastic systems. Statist. Inf. Stochast. Proces. 3, 173-182. 

[70] M. L. Kleptsyna and A. Le Breton (2002), Statistical analysis of the fractional Ornstein-Uhlenbecktype process, 
Statistical Inference for Stochastic Processes 5(3), 229-242. 

[71] F. Klebaner (2005), Introduction to Stochastic Calculus with Application, Imperial College Press, Second Edi- 
tion, 2005. 

[72] S. Kogon, D. Williams (1998), Characteristic function based estimation of stable distribution parameters, in A 
practical guide to heavy tails, M. T. R. Adler R. Feldman, Ed. Berlin: Birkhauser, 31 1-335. 

[73] A.N. Kolmogorov (1940), Wienersche spiralen und einige andere interessante kurven im hilbertschen raum, 
C.R.(Doklady) Acad. URSS (N.S) 26, 115-118, 1940. 

[74] I. A. Koutrouvelis (1980), Regression-type estimation of the parameters of stable laws. Journal of the American 
Statistical Association 75, 918-928. 



26 JIARUI YANG AND JINQIAO DUAN 

[75] H. Kunita (2004), Stochastic differential equations based on Levy processes and stochastic flows of diffeomor- 
phisms, Real and stochastic analysis (Eds. M. M. Rao), 305-373, Birkhuser, Boston, MA. 

[76] E. E. Kuruoglu (2001), Density parameter estimationof skewed a-stable distributions, Singnal Processing, IEEE 
Transactions on 2001, 49(10): 2192-2201. 

[77] Yu. A. Kutoyants (1984), Parameter estimation for diffusion type processes of observations, Statistics 15(4), 
541-551. 

[78] A. Le Breton (1998), Filtering and parameter estimation in a simple linear model driven by a fractional Brown- 
ian motion, Stat. Probab. Lett. 38(3), 263-274. 

[79] A. Le Breton (1976), On continuous and discrete sampling for parameter estimation in diffusion type processes, 
Mathematical Programming Study 5, 124-144. 

[80] R. S. Lipster and A. N. Shiryaev (1977), Statistics of Random Processes, Springer, New York, 1977. 

[81] X. Liu, J. Duan, J. Liu and P. E. Kloeden (2009), Synchronization of systems of Marcus canonical equations 
driven by a-stable noises, Nonlinear Analysis - Real World Applications, to appear, 2009. 

[82] A. W. Lo (1991), Long-term memory in stock market prices, Econometrica 59, 1279-1313. 

[83] B. B. Mandelbrot and J. R. Wallis (1969), Computer experiments with fractional Gaussian noises, Water Re- 
sources Research 5, 228-267. 

[84] B. B. Mandelbrot and J. W. Van Ness (1968), Fractional Brownian motions, fractional noises and applications, 
SI AM Rev. 10, 422-437. 

[85] X. Mao (1995), Stochastic Differential Equations and Applications, Horwood Publishing, Chichester. 

[86] B. Maslowski and B. Schmalfuss (2005), Random dynamical systems and stationary solutions of differential 
equationsdriven by the fractional Brownian motion, Stoch. Anal. Appl. 22(6), 1577 - 1607. 

[87] J. H. McCulloch (1986), Simple consistent estimators of stable distributions, Communications in Statistics- 
Simulation and Computation 15, 1109-1136. 

[88] Y. S. Mishura (2008), Stochastic Calculus for Fractional Brownian Motion and Related Processes, Springer, 
Berlin. 

[89] F. Moss and P. V. E. McClintock (eds.), Noise in Nonlinear Dynamical Systems. Volume 1: Theory of Contin- 
uous Fokker-Planck Systems (2007); Volume 2: Theory of Noise Induced Processes in Special Applications 
(2009); Volume 3: Experiments and Simulations (2009). Cambridge University Press. 

[90] J. P. Nolan (2007), Stable Distributions - Models for Heavy Tailed Data, Birkhause, Boston, 2007. 

[91] I. Nourdin and T. Simon (2006), On the absolute continuity of Levy processes with drift, Ann. Prob. 34(3), 
1035-1051. 

[92] I. Norros, E. Valkeila and J. Virtamo (1999), An elementary approach to a Girsanov formula and other analytical 
results on fractional Brownian motions, Bernoulli 5(4), 571-587. 

[93] D. Nualart (2003), Stochastic calculus with respect to the fractional Brownian motion and applications, Con- 
temporary Mathematics, 336, 3-39. 

[94] B. Oksendal (2005), Applied Stochastic Control Of Jump Diffusions, Springer- Verlag, New York. 

[95] B. Oksendal (2003), Stochastic Differenntial Equations, Sixth Ed., Springer- Verlag, New York. 

[96] B. Oksendal, F. Biagini, T. Zhang and Y. Hu (2008), Stochastic Calculus for Fractional Brownian Motion and 
Applications. Springer. 

[97] T. N. Palmer, G. J. Shutts, R. Hagedorn, F. J. Doblas-Reyes, T. Jung and M. Leutbecher (2005), Representing 
model uncertainty in weather and climate prediction, Annu. Rev. Earth Planet. Sci. 33, 163-193. 

[98] A. Papoulis (1984), Probability, Random Variables, and Stochastic Processes, McGraw-Hill Companies, 2nd 
edition. 

[99] N. D. Pearson and T. Sun (1994), Exploiting the conditional density in estimating the term structure: an appli- 
cation to the Cox, Ingersoll and Ross model, The Journal of Finance 49(4), 1279-1304. 

[100] A. R. Pedersen (1995), Consistency and asymptotic normality of an approximate maximum likelihood estimator 
for discretely observed diffusion processes, Bernoulli 1(3), 257-279. 

[101] C. K. Peng, V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, and A. L. Goldberger(1994), Mosaic organization 
of DNA nucleotides. Phys. Rev. E 49, 1685-1689. 

[102] S. Peszat and J. Zabczyk (2007), Stochastic Partial Differential Equations with Levy Processes, Cambridge 
University Press, Cambridge, UK. 



QUANTIFYING UNCERTAINTIES IN COMPLEX SYSTEMS 



27 



[103] B.L.S. Prakasa Rao (1999), Statistical Inference for Diffusion Type Processes, Arnold, London. 

[104] B.L.S. Prakasa Rao (1999), Semimartingales and their Statistical Inference, Chapman & Hall/CRC. 

[105] B. L. S. Prakasa Rao (2003), Parametric estimation for linear stochastic differential equations driven by frac- 



tional Brownian motion, http://www.isid.ac.in/~statmath/eprints 



[106] S. Press (1972), Estimation of univariate and multivariate stable distributions, Journal of 'the Americal Statistical 
Association 67, 842-846. 

[107] P. E. Protter (2005), Stochastic Integration and Differential Equations, Springer- Verlag, New York, Second 
Edition. 

[108] B. L. Rozovskii (1990), Stochastic Evolution Equations, Kluwer Academic Publishers, Boston. 
[109] P. M. Robinson (1977), Estimation of a time series model from unequally spaced data, Stock. Proc. Appl. 6, 
9-24. 

[110] G. Samorodnitsky, M. S. Taqqu (2008), Stable Non-Gaussian Random Processes- Stochastic Models with Infi- 
nite Variance. Chapman & Hall/CRC. 

[Ill] K. Sato (1999), Levy Processes and Infinitely Divisible Distrributions, Cambridge University Press, Cambridge, 
UK, 1999 

[112] H. Scher, M. F. Shlesinger and J. T. Bendler (1991), Time-scale invariance in transport and relaxation, Phys. 
Today 44(1), 26-34. 

[113] D. Schertzer, M. Larcheveque, J. Duan, V. Yanovsky and S. Lovejoy (2000), Fractional Fokker-Planck equation 
for non-linear stochastic differential equations driven by non-Gaussian Levy stable noises, J. Math. Phys. 42, 
200-212. 

[114] M. F. Shlesinger, G. M. Zaslavsky and U. Frisch (1995), Levy Flights and Related Topics in Physics, Lecture 

Notes in Physics, Springer- Verlag, Berlin. 
[115] M. Sorensen (1999), On asymptotics of estimating functions, Brazillian Journal of Probability and Statistics 

13, 111-136. 

[116] D. W. Stroock and S. R. S. Varadhan (1979), Multidimensional Diffusion Processes, Springer Verlag, Berlin. 

[117] T. H. Solomon, E. R. Weeks, and H. L. Swinney (1993), Observation of anomalous diffusion and Levy flights 
in a two-dimensional rotating flow, Phys. Rev. Lett. 71(24), 3975 - 3978. 

[118] M. S. Taqqu, V. Teverovsky, and W. Willinger (1995), Estimators for long-range dependence: an empirical 
study, Fractals, 3(4), 785-798. 

[119] G. A. Tsihrintzis and C. L. Nikias (1995), Fast estimation of the parameters of alpha-stable impulsive interfer- 
ence using asymptotic extreme value theory, ICASSP-95, 3, 1840-1843. 

[120] N. G. Van Kampen (1987), How do stochastic processes enter into physics? Lecture Note in Mathe. 1250/1987, 
128-137. 

[121] N. G. Van Kampen (1981), Stochastic Processes in Physics and Chemistry, North-Holland, New York. 
[122] E. Waymire and J. Duan (Eds.) (2005), Probability and Partial Differential Equations in Modern Applied Math- 
ematics, Springer- Verlag. 

[123] D. S. Wilks (2005), Effects of stochastic parameterizations in the Lorenz '96 system, Q. J. R. Meteorol. Soc. 
131, 389-407. 

[124] P. D. Williams (2005), Modeling climate change: the role of unresolved processes, Phil. Trans. R. Soc. A 363, 
2931-2946. 

[125] E. Wong and B. Hajek (1985), Stochastic Processes in Engineering Systems, Spring-Verlag, New York. 

[126] W. A. Woyczynski (2001), Levy processes in the physical sciences, In Levy processes: theory and applications, 

O. E. Barndorff-Nielsen, T. Mikosch and S. I. Resnick (Eds.), 241-266, Birkhauser, Boston, 2001. 
[127] A. M. Yaglom (1958), Correlation theory of processes with random stationary nth increments, AMS Transl. 

2(8), 87-141. 

[128] Z. Yang and J. Duan (2008), An intermediate regime for exit phenomena driven by non-Gaussian Levy noises, 

Stochastics and Dynamics 8(3), 583-591. 
[129] F. Yonezawa (1996), Introduction to focused session on 'anomalous relaxation, J. Non-Cryst. Solids 198-200, 

503-506. 

[130] N. Yoshida (2004), Estimation for diffusion processes from discrete observations, J. Multivariate Anal. 41(2), 
220-242. 



28 JIARUI YANG AND JINQIAO DUAN 

(Jiarui Yang and Jinqiao Duan) Department of Applied Mathematics, Illinois Institute of Technology, Chicago, IL 
60616, USA 

E-mail address: J . Yang: jyang31@iit.edu; J. Duan: duan@iit.edu 



i r 



200 300 400 500 600 700 800 900 

N 



